Molecular Systems Biology

Interpreting the functional consequences of genetic variation, especially rare missense variants, remains a significant challenge in human genetics. Computational variant effect predictors (VEPs) and multiplexed assays of variant effects (MAVEs) provide complementary approaches, with VEPs offering scalable predictions and MAVEs delivering detailed empirical measurements. However, MAVEs are resource intensive and cannot yet be applied broadly across the proteome, making it important to identify proteins where experimental mapping will be most informative. We hypothesised that MAVEs should be particularly valuable for proteins where computational predictors disagree, as such disagreement may highlight mechanistic blind spots. To test this, we analysed predictions from ten distinct VEPs across more than 13,000 human proteins and quantified inter-predictor concordance. We observed substantial variability across proteins in the degree of agreement across predictors and investigated structural, functional and gene-level features associated with this variation. We find that inter-VEP concordance showed no relationship with agreement to experimental MAVE data. If predictor agreement reflected how intrinsically predictable a protein is, these quantities would be expected to correlate. Their decoupling instead suggests that MAVEs may provide orthogonal information to VEPs, supporting the use of inter-VEP disagreement to prioritise proteins where experimental data will be most informative. We therefore propose using inter-VEP disagreement as a practical strategy to prioritise proteins for experimental characterization. Focusing on proteins with low predictor concordance should maximise the informational value of new MAVEs, and improve variant interpretation in both research and clinical contexts.

The ability to understand and affect the course of complex, multi-system diseases like diabetes has been limited by a lack of well-designed, high-quality and large multimodal datasets. The NIH Bridge2AI AI-READI project (aireadi.org) aims to address this shortfall by generating an AI-ready dataset to support AI discoveries in type 2 diabetes mellitus (T2DM). This manual of procedures provides a detailed description of the AI-READI protocol.

Functional health is centered on five domains of Intrinsic Capacity (IC): locomotion, cognition, vitality, psychological and sensory capacity. Therefore, measuring IC at the domain-specific level is the cornerstone for developing preventive interventions to help individuals preserve their independence. In this study, we used 63 clinical features from the UK Biobank to develop IC age, an 18-year mortality risk estimator that approximates an individuals biological age associated with the decline of each IC domain. By establishing proteomic surrogates of IC age, we find immune system activation across domains and provide a proteomic framework that may facilitate scalable monitoring of functional health decline.

PDAC exhibits rapid chemoresistance, yet how drug-tolerant states arise remains unclear. Existing approaches miss how network topology evolves across cell-state transitions under drug pressure. A 3D PANC-1 tissue model on decellularized intestinal matrix was used for scRNA-seq across four conditions (control, GEM, TGF-{beta}1, GEM+TGF-{beta}1). Pseudotime trajectory inference was combined with dynamic PPI network analysis. Findings were cross-examined in a PDAC atlas (726,107 cells, 231 patients; Loveless et al., 2025). GEM resistance involved E2F1, mTOR, CDK1, AURKA, TPX2, TOP2A, and BIRC5. TGF-{beta}1 drove EMT resistance via KRAS, glycolysis, and hypoxia, inducing SPOCK1, MBOAT2, COL5A1, ADAMTS6, THBS1, and FN1. Trajectory-coupled network analysis revealed an emergent bottleneck when G1[->]S and TGF-{beta}1-induced EMT co-occurred: CDK1 centrality spiked selectively, with CDKN1A as critical regulator. This CDK1-CDKN1A-WEE1 axis defines an "S-phase persistence" state enriched for GEM survivors. Atlas cross-examination confirmed 8.7-fold metastatic enrichment of triple-positive cells and EMT-cell-cycle coupling. Trajectory-coupled network topology analysis identifies CDK1-CDKN1A-WEE1 as a chemoresistance bottleneck corroborated in 726,107 patient cells. The framework generalizes to drug resistance across cancer types.

Wastewater (WW) genomics can track SARS-CoV-2 circulation beyond clinical testing, but its ability to reflect clinical diversity and capture severity-linked mutations remains unclear. Here, we integrated 845 clinical genomes and 22 wastewater genomes from Valencia, Spain, across matched metropolitan and hospital catchments. We compared matched WW and clinical sequencing for lineage and mutation surveillance at two levels: metropolitan and hospital. Then, we tested WW sensitivity to detect mutations statistically associated with hospitalization status in regional (n = 4,843), national (n = 10,052) and supranational (n = 39,099) clinical datasets. WW surveillance captured the dominant Omicron background when collapsing lineages into parental lineages constellation but had limited sensitivity for fine-scale sublineage diversity. Performance was strongly catchment dependent: metropolitan wastewater best represented broader community circulation, whereas hospital wastewater was noisier but detected KP.3 months before its appearance in routine metropolitan clinical surveillance. Across clinical datasets, hospitalisation-associated substitutions showed limited reproducibility, although the national and supranational analyses converged on receptor-binding-domain substitutions D405N, K417N and R408S. Networks showed coupling between G252V in NTD with those RBD substitutions involved in immune escape and receptor engagement. Finally, integrating regional to supranational GWAS with interaction networks and wastewater detection prioritised mutations supported by at least two independent association layers, that includes mutations in the Spike, especially in RBD, and the wastewater-exclusive candidate S:V445P, which was missed by contemporaneous clinical sequencing. Overall, WW genomics preferentially recovers the common mutational backbone of SARS-CoV-2 circulation and can highlight important changes missed by clinical sampling, making it a complementary tool for real-time prioritisation of viral evolutionary change.We found partial overlap in lineage composition between WW and clinical samples, with higher overlap at the metropolitan (50%), than at the hospital level (30%). Conversely, we found a slightly higher overlap of individual mutations between WW and clinical samples at the hospital level (20%) than at the metropolitan area (16%), but shared mutations in both datasets were enriched in the Spike gene. Clade composition did not differ between 216 hospitalised and 528 non-hospitalised cases at regional level. Using GWAS and Hierarchical Lasso analysis, we detected mutations associated with hospitalization status in three different datasets: regional, national and worldwide, with little overlap between them. Although few variants replicated across cohorts, the overlap between the Spain and global analyses was statistically enriched and centred on RBD substitutions (D405N, K417N, R408S). Multiple integration of genomic association results prioritised 34/191 wastewater mutations (16 in Spike), including one mutation only detected in wastewater missed by routine clinical surveillance. Wastewater sequencing tracked dominant Omicron waves but performance varied by catchment; integrating clinical association results with interaction network modelling helped prioritise and interpret wastewater-detected mutations.

Type 2 diabetes case reports describe complex clinical courses, but their timelines are often expressed in language that is difficult to reuse in longitudinal modeling. To address this gap, we developed a textual time-series corpus of 136 PubMed Open Access single-patient case reports involving glucagon-like peptide 1 receptor agonists, with clinical events associated with their most probable reference times. We evaluated automated LLM timeline extraction against gold-standard timelines annotated by clinical domain experts, assessing how well systems recovered clinical events and their timings. The best-performing LLM produced high event coverage (GPT5; 0.871) and reliable temporal sequencing across symptoms (GPT5; 0.843), diagnoses, treatments, laboratory tests, and outcomes. As a downstream demonstration, time-to-event analyses in diabetes suggested lower risk of respiratory sequelae among GLP-1 users versus non-users (HR=0.259, p<0.05), consistent with prior reports of improved respiratory outcomes. Temporal annotations and code will be released upon acceptance.

Predicting biological responses to perturbations such as stress, nutrition, or pharmaceuticals could transform healthcare and biosciences. The task is challenging because complex interactions among many factors interactively drive biological variation, and thereby alter responses. However, metabolism integrates these drivers of variation, suggesting that emergent biological states may be reflected in aggregate chemical states. We define such states as chemotypes and test their predictive utility. Using Fourier-transform infrared (FTIR) spectroscopy coupled with machine learning in Drosophila melanogaster, we show that chemotypes serve as proxies for biological variation driven by sex, genotype, nutrition and age, and - critically - predict among-population variation in stress response. These findings indicate that chemotypes provide a computable and integrative representation of organismal biology, predicting genotype, phenotype, and response to perturbation.

The spatial arrangement of cells within a tissue microenvironment shapes their interactions and cell states, which are essential for tissue development, homeostasis, and disease. Spatial -omics technologies can precisely map the location of each cell within complex tissue structures, while also profiling their protein content and transcriptional diversity. Various approaches have been developed to analyze spatial patterns of cell aggregation, repulsion, or random distribution within tissues. However, differences in cell morphology within a tissue can introduce significant bias. Cell size in particular is not accounted for and introduces challenges when quantifying the aggregation of cells or their molecular features. To overcome such limitations, we present ClumPyCells: a statistical framework that measures cell and marker aggregation within tissue while correcting for size morphology. ClumPyCells enables interpretation of cell aggregation, bypassing interfering cell types or tissue regions unrelated to the desired spatial correlation. We demonstrate the capabilities of ClumPyCells across several tumor types, including melanoma and colorectal cancer, and spatial -omics technologies such as spatial transcriptomics and proteomics, while benchmarking how cell-size differences contribute to misinterpretations. By correcting for disruptive cell types within a region of interest, ClumPyCells will determine new tissue patterns and structures without morphological interference.

The gut modulates systemic health, influencing immune, neurological, and cardiovascular processes. While fecal sequencing of microbial nucleic acids provides a non-invasive view of microbial composition, sensitive measurement of host-derived signals in stool remains limited. Here we introduce DIGEST (Digital Immunoassay for Gut-Environment Single-molecule Targets), an ultrasensitive digital immunoassay that quantifies proteins in fecal extracts to attomolar levels. In mice, longitudinal profiling during a high-fat diet perturbation revealed coordinated host cytokine responses that occurred within 24 hours, with sustained elevation after diet withdrawal, enabling non-invasive tracking of within-subject immune dynamics. Application of DIGEST to quantify a panel of host inflammatory cytokines in patients with inflammatory bowel disease distinguished active ulcerative colitis from quiescent disease and non-IBD controls (AUC=0.98). In advanced melanoma patients receiving PD-1 blockade, pretreatment fecal IL-23 concentrations discriminated responders from non-responders with an AUC of 0.87. Together, these results establish DIGEST as a generalizable platform for sensitive, non-invasive quantification of host protein activity at the gut interface, with broad applications in basic science discovery, disease surveillance, and therapy response prediction.

Lipoprotein composition is altered in sepsis, and supplementation with high-density lipoproteins has been reported to improve outcomes in experimental settings. In this study, we aimed to investigate the nature and inter-individual variability in the lipoprotein proteome to inform risk stratification and opportunities for precision medicine approaches. In a large proteomic dataset including 1134 patients (1781 samples) with sepsis and 149 healthy volunteers, we analysed 18 protein components of lipoproteins. We characterise heterogeneity of the lipoprotein proteome, defining three step-wise sub-phenotypes associated with increasing disease severity, one close to health, then an early phase patient group showing increased abundance of proteins that integrate HDL under inflammatory conditions (SAA1 and SAA2), then a group with decreased abundance of proteins that are components of HDL under healthy conditions that was associated with higher organ failure intensity (SOFA score) and increased mortality. We developed and externally validated a quantitative score reflective of lipoproteins alterations in sepsis, and machine learning predictive models to predict the LP class, advancing future individualised lipoproteins-based therapeutics in sepsis.

Tissues are organized with interacting multicellular organizational units whose interfaces and transitions shape function in health and disease. Current spatial-omics analyses typically assign cells to a single cellular neighborhood--ignoring natural gradients, heterogeneity, and borders. Here we present MINGL (Mixture-based Identification of Neighborhood Gradients with Likelihood estimates), a probabilistic framework that converts existing neighborhood annotations into continuous measures of tissue architecture. MINGL models each cell by multi-membership probabilities across hierarchical organizational units and uses these probabilities to identify enriched cells at interfaces between units, constructs interaction networks across hierarchical scales, quantifies compositional gradient transitions, measures context-specific composition heterogeneity, and provides a starting point for neighborhood resolution. Across multiple spatial-omic datasets spanning melanoma, healthy intestine, and Barretts Esophagus progression, MINGL detected innate immune-enriched interfaces at tumor and anatomical interfaces, plasma cell niches linking cellular neighborhoods, distinct regimes of sharp and gradual transitions between organizational states, and disease-associated neighborhood remodeling. By treating neighborhood assignment uncertainty as a biological signal rather than noise, MINGL unifies discrete and continuous representations of tissue organization and makes tissue architecture measurable, comparable, and scalable across biological scales and spatial-omics platforms.

The RAS/RAF/MEK/ERK pathway has been extensively studied for its roles in physiology and disease, yet a systems-level view of its interplay with the broader cellular context is lacking, limiting insight into paralog-specific roles, mutation-driven rewiring, and cross-pathway integration. To address this, we simultaneously mapped the interactomes of all pathway components and their activating and inactivating mutations, identifying 2,500 high-confidence interactors, 88% of which previously unreported, and generating a high-confidence, comprehensive reference map of the pathway in resting and activated states. The pathway behaves like a "molecular organism" that connects to other signaling complexes, operates from distinct subcellular sites, and integrates into networks relevant to physiology and disease. Dataset exploration uncovered paralog- and mutation-specific interaction signatures and pathway-level crosstalk, including previously unknown connections to mRNA metabolism and WNT signaling. The interactomes available at RAStoERK.univie.ac.at provide a foundational community resource for unbiased discovery, targeted mechanistic studies, and potential therapeutic exploration. HighlightsO_LIHigh-resolution interactome map of the RAS to ERK pathway in resting and active states C_LIO_LIState-dependent pathway-level interactions with signaling and functional modules C_LIO_LIMutation-driven interactome rewiring and integration into pathophysiological networks C_LIO_LIA resource to discover pathway crosstalk, disease links, and actionable interventions C_LI

Our bodies have evolved to maintain homeostasis through regulatory systems that continuously adapt to keep physiological processes within a normal range. From this perspective, complex chronic disease can be understood as a breakdown of compensatory mechanisms, resulting in loss of homeostasis. Here we propose that adaptive receptor expression dynamics may serve as one such compensatory mechanism, increasing receptor surface expression when external ligand is insufficient, and clearing it when signaling is excessive. To explore this, we adapt a previously published agent-based model and simulate it under a range of scenarios. We find that the system of adaptive receptor expression is robust to oscillatory perturbations but not to chronic stress. We propose that receptor turnover dynamics may be better understood as an adaptive, environmentally responsive process rather than a fixed biological property, and that in some cases, disease manifests only after compensatory mechanisms have been pushed past their limits. We conclude with a discussion of implications for understanding complex chronic diseases, for thinking about epigenetic and mutational change as escalating layers of adaptation, and for how we model receptor dynamics in the context of receptor-mediated drug activity.

Antibiotic resistant pathogens are an increasing public health threat, as development of novel therapeutics is outpaced by resistance emergence and dissemination. Approaches to slow down or even revert antibiotic resistance are necessary to maintain efficacy of both existing and new antibiotics. Such approaches exploit the fitness cost of resistance elements, but have largely relied on assessing this cost in laboratory conditions that poorly reflect the native context in which pathogens reside. Here we present a method that allows to investigate the influence of personalized human gut microbiota compositions on the competitive fitness of antibiotic resistant pathogens. Using fecal matter-derived microbiomes we identify a specific community that selects for a carbapenem-resistant Klebsiella pneumoniae strain. This selective advantage is due to mutations arising in a LacI-type transcriptional regulator, GlyR. We show that upregulation of the downstream glycoporin GlyP is causing the effect. By deconvoluting the microbiome composition, we identify a focal E. coli strain as a central driver of the selection, which is further modulated by other microbiota members. We demonstrate that the selective advantage is due to carbohydrate competition, and in particular for glycerol-containing compounds. Importantly, glyR mutations are under strong positive but conditional selection in clinical K. pneumoniae isolates. This implies a reduced competitiveness in other environments, which we experimentally validate in vitro. Overall, this study offers a path to identify microbiome-specific interactions that modulate the competitiveness of antibiotic resistant pathogens.

Pancreatic beta cells have the unique function of synthesizing and secreting high amounts of the inhibitory neurotransmitter {gamma}-aminobutyric acid (GABA). The mechanism of GABA secretion, whether vesicular or channel-mediated, is debated. Our study reveals surprising temporal complexity in the pattern of islet GABA secretion. We used insulin secretion modulators to demonstrate that GABA release is not directly correlated with insulin secretion. VGAT reporter mice also showed that beta cells do not express the requisite vesicular GABA transporter (VGAT) for vesicular GABA release. Instead, GABA is secreted from the cytosol in pulses by the LRRC8A/D isoform of the volume regulatory anion channel (VRAC). We further demonstrate the dynamic coordination of GABA release with calcium influx in beta cells and dependence on beta cell depolarization. These results suggest a model where GABA is released during the peaks of beta cell calcium oscillations to provide feedback which strengthens and reinforces the oscillation waveform.

Synonymous codons encode the same amino acid but can differ in their usage and translational properties. In previous work we reported statistical differences in backbone dihedral angle distributions associated with synonymous codons in the Escherichia coli proteome. This finding has been questioned due to concerns regarding the statistical methodology used. Here we revisit the dataset using corrected statistical procedures and alternative statistical tests. Across multiple frameworks, the real dataset consistently shows an excess of small p-values relative to randomized controls, indicating detectable codon-associated differences in backbone conformation.

Antimony is a human-readable language for defining and sharing models developed by the systems biology community. It enables scientists to describe biochemical networks with a simple syntax, while supporting seamless conversion to and from the Systems Biology Markup Language (SBML) community standard. Since Antimonys original release, both SBML and modeling practices have evolved significantly, creating a need to update Antimony to maintain its standards compliance and practical relevance. In this paper, we introduce Antimony 3, a comprehensive update that formalizes its cumulative improvements and extends its support for SBML Level 3 Core and Flux Balance Constraints (FBC), Distributions, Layout, and Render packages. Antimony 3 enables model specifications that combine kinetic reactions with flux balance analysis, represent uncertainty using probability distributions, add biological context through annotations, and define publication-ready visualizations, all within a unified plain-text format. Antimony 3 is delivered as a lightweight C/C++ library with a stable C API. It is available through official bindings for Python, Julia, and JavaScript/WebAssembly, as well as a cross-platform desktop GUI, which enables straightforward use across scripting environments, desktop applications, and browser-based tools. Antimony 3 is released as open-source software under the BSD 3-Clause License and is available at https://github.com/sys-bio/antimony. Author SummaryBiological models are typically stored in standardized formats that ensure compatibility across different software tools, but these formats rely on verbose, machine-readable syntax that is difficult for humans to write or inspect directly. Antimony addresses this challenge by providing an intuitive, text-based language for defining biological models that can be automatically converted to and from the Systems Biology Markup Language (SBML). Since Antimonys original release in 2009, the SBML standard and common modeling workflows have expanded significantly. We developed Antimony 3 to support these advances, enabling researchers to write a single human-readable text file that defines reaction networks, constraint-based objectives, uncertainty in parameters and initial conditions, semantic annotations linking to biological databases, and model diagrams. Antimony 3 is provided as open-source software with broad support across computational environments, making it accessible to researchers in a wide range of workflows.

Spatiotemporal organisation of biological molecules is a key driver of cellular processes, including many post-transcriptional epigenetic processes. The germline-specific germ granules are biomolecular condensates that act as hubs for mRNA and small RNA processing and are core regulators of germline gene expression programming. Germ granules have been studied extensively in C. elegans, and recent developments have led to many subdivisions of the germ granule into specialised compartments. Rapid advancements in microscopy and protein-protein interaction (PPI) screening techniques have produced a large amount of data towards characterising the localisation of proteins to specific granules. However, common methods used to probe PPIs are limited in their ability to robustly detect valid interactions, especially the multivalent and sometimes transient ones observed in granule environments. Here we perform a meta-analysis of granule protein interaction screens. While these experiments generally enrich for proteins matching the profile of granule-associated proteins, we find that when considering screens individually, reproducibility is surprisingly low, highlighting not only the variability inherent in these methods but also the dynamic nature of the PPI networks present in granules. We developed an algorithm to provide a measure of each proteins association with specific granules across various experiments. By further clustering and investigation of the resulting score matrix, we demonstrate the power of this holistic approach to provide deeper insights into germ granule organisation and highlight novel can provide a resource to better inform future investigations into granules and their constituent proteins.

Tumor-induced bone disease (TIBD) arises from a complex interplay between metastatic cancer cells and the bone microenvironment, creating a self-reinforcing "vicious cycle" of bone destruction and tumor growth. Experimental evidence from our group (Buenrostro et al., Bone 113:77-88, 2018) suggests that tumor cells in the bone microenvironment early in disease rely more heavily on bone-derived growth factors, such as transforming growth factor-{beta} (TGF-{beta}), to sustain proliferation than tumor cells late in disease, which may grow independently of these factors. Here, we integrate a mechanistic, population-dynamics model of tumor-bone interactions with in vivo data to test the hypothesis that inhibiting bone resorption suppresses growth of non-adapted but not bone-adapted tumors. The model includes key regulators of TIBD, including TGF-{beta}-driven tumor proliferation, parathyroid hormone-related protein (PTHrP) secretion, and osteoblast (OB)-osteoclast (OC) coupling. Parameter calibration using data from mice injected intratibially with parental (non-adapted) and bone-adapted breast cancer cells reveals distinct parameter values for each tumor type. Bone-adapted cells exhibit a higher basal division rate and reduced sensitivity to TGF-{beta}-mediated stimulation, whereas parental-derived tumor cells depend more strongly on TGF-{beta} and secrete PTHrP at higher rates to compensate for their slower growth. Model simulations reproduce the greater bone loss observed experimentally for bone-adapted tumors and predict that, for non-adapted tumors, bone destruction results from a slower but meaningful rise in OC activity and a possible moderate decline in OBs. Simulated treatment of bone-adapted tumors with the bisphosphonate zoledronic acid stabilizes bone density but has limited or highly variable effects on tumor growth. These results suggest that OC inhibition alone may be insufficient to restrain tumor expansion once tumors have adapted to the bone microenvironment. Together, these findings support the hypothesis that tumor adaptation to the bone microenvironment governs dependence on bone-derived growth factors and response to OC-targeted therapy, underscoring the value of mechanistic modeling for elucidating tumor-bone interactions and guiding tumor-type-specific treatment strategies for TIBD.

SummaryExtracellular vesicles (EVs) are bilayer vesicles that carry a diverse cargo of molecules, such as nucleic acids, proteins and metabolites. These EVs can be transported throughout the organism to specific recipient tissues. For this reason, EVs have been recognized as pivotal mediators of cell-to-cell communication (CCC). Importantly, alterations in EV-mediated communication have been linked to pathological processes, further highlighting their biological relevance. However, the in silico exploration of the functional effects of EV cargo in recipient tissues remains limited due to the lack of dedicated tools that can be applied to EV omics datasets. Most current bioinformatics tools for assessing CCC rely on ligand-mediated communication and therefore cannot be used to explore EV-mediated communication. To address this gap, we developed EV-Net, a bioinformatics tool designed to explore the effects of EV cargo on recipient tissues. EV-Net was built by adapting NicheNet, a CCC bioinformatics tool that relies on ligand-receptor mediated communication, for the analysis of EVs proteomics and RNA-seq data. The EV-Net framework enables the identification and prioritization of EV cargo molecules with high regulatory potential in a recipient tissue of interest. This prioritization facilitates the systematic translation of EV cargo profiles into testable biological hypotheses. Availability and documentationThe source code of EV-Net is stored in GitHub https://github.com/torrejoNia/EV-Net alongside instructions on how to install it. Comprehensive tutorials and additional documentation are available at https://torrejonia.github.io/EV-Net/. The datasets used in the use cases were deposited in Zenodo. The corresponding Zenodo links are provided in the tutorials for each use case. This software is distributed under a GLP3 licence.